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Abstract 

In this paper we prove the existence of the high-energy tails for electron distri- 
bution function of the Boltzmann equation for semiconductors, in the stationary 
and homogeneous regime, in the analytic band approximation and scattering 
with acoustic and optical phonons and impurities. We also prove numerically 
that the tail is a global maxwellian in the parabolic band approximation, and 
in the quasi parabolic band case, a power law of the global maxwellian. 
keywords Boltzmann-Poisson system for semiconductors, Monte Carlo Method, 
Semiconductors. 
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1 Introduction 

In recent years, the dimensions of MOS devices have been scaled down reaching sub- 
micrometric order. As a result, large electric fields near the drain region generate hot 
or energetic electrons. Hot electrons can heat the device with important consequences 
for long-term reliability, or can be injected into the oxide creating an instability in 
the device performance. Moreover, injection of hot electron onto a floating gate has 
been widely used as the writing mechanism in EPROM and flash EEPROM device. 
To model these effects is a difficult task, because they depend upon the small fraction 
of carriers with energy above certain thresholds and hence the knowledge of the tail 
of the electron energy distribution (hereafter EED) is required. 
The natural framework for describing these regimes is the Boltzmann Transport Equa- 
tion (BTE) with the adequate scattering mechanisms, coupled with the Poisson equa- 
tion. To solve the BTE is an hard task: recently deterministic solvers of the BTE, 
based of finite difference schemes, have been introduced [IJ but their computational 
efforts is comparable to that of the Direct Simulation Monte Carlo. 



Another alterative is to use Extended hydrodynamic models, based on the higher 
order moments of the BTE j2] , but their sensitivity to the tail of the EED is an open 
question. 

At the moment, the Monte Carlo (MC) solution of the BTE represents the most 
accurate approach for describing the tail of the EED, because the correct scattering 
mechanisms as well as the bands structure is taken into account. Unfortunately, 
the MC approach is CPU consuming, hence possibly unsuitable for extensive use 
in technology development. What we know from MC simulations is that the EED 
shape, for high energies, is maxwellian [HJ Hj (for parabolic bands and acoustic and 
optical phonons scattering), whereas the electron-electron interaction modifies this 
asymptotic behaviour enhancing the tail O E] . 

The understanding of these phenomena has been the subject of a wide debate: ac- 
cording to some authors such tails reflect the distribution at the device boundaries 
due to ballistic transport j3IH], or has been caused by 'superlucky' electrons which 
have absorbed more phonons than they have emitted [SI. 

A lot of effort has also been put into the development of analytical expression for 
the whole EED. Heuristic models for the EED are based on expansions around a 
global maxwellian [^j , or polynomial in the exponential function [21021! or a combi- 
nation of a global maxwellian and a power law of a global maxwellian (also known as 
'stretched exponential') [TT] in order to describe 'cold' and 'hot' electron populations. 
The parameters appearing in these EEDs can be determined by considering some 
approximation of the BTE and by fitting the resulting model to MC data. 
The problem of the description of the high-energy tails has been tackled in a different 
context as in kinetic models for granular flows. A variety of analytical techniques have 
been developed such as by comparing the 'gain' and 'loss' terms in the BTE ^21, or 
by means of a rigorous pointwise lower estimate |T3] , or by means of weighted tail 
control Jl] , or by using the so-called inelastic Maxwell models [T^l ITE] . 
The aim of this paper is to tackle the high-energy tail problem rigorously in the BTE 
framework, by using an analytical technique adapted from granular media which is 
based on moments estimate In this way, we are able to prove the existence of 
the EED tail which has the form of a 'stretched exponential'. The paper is organized 
as follows. In Section 2 the basic equations are written. In Section 3 the 'stretched 
exponential' for the tail is introduced and we formulate the main results. Section 4 
is devoted to the moment inequalities. In Section 5 we present the proof of the main 
theorem. In Section 6 simulation results obtained by means of Direct Simulation 
Monte Carlo are shown and conclusion are drown in Section 7. 



2 Basic Equations 

The Boltzmann transport equation (BTE) for one conduction band ^3 QHl is 

|^ + v(k)-Vx/-|E-VJ=Q[/] , (1) 

where the unknown function f{t,x, k) represents the probability density of finding an 
electron at time t in the position x = {xi,X2, x^) with the wave- vector k = {ki, k2, ks), 
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and q is the absolute value of the electron charge. The domain of k is the first Brillouin 
zone B. In the neighborhood of the band minimum the dispersion relation can be 
considered approximately quasi-parabolic: 

e(k)[l + ae(k)]=7 = ^, (2) 

where m* denotes the effective electron mass, which is 0.32 rrie (free electron mass) in 
silicon, a the non parabolicity factor (0.5 eV^~^),and h the Planck constant divided 
by 27r. With this approximation the Brillouin zone B coincides with M^. 
From eq.(j2I), the electron group velocity v = (fi,f2,t'3) writes 

/, N 1 ^ ^k 1 
V k = -Vk^ = , . 3 

The parabolic band approximation is obtained from the previous equations with a=0. 
The electric field E(t,x) = {Ei, E2, E3) satisfies the Poisson equation 

A{e<P) = q[n{t,x) - Nd{x) + N^ix)] 
E = -Vx</> 

where 0(t, x) is the electric potential, A''^ and Na, respectively, are the donor and 
acceptor densities, e the dielectric constant, n the electron density given by 



n(t,x)= / /(t,x,k). 
In general, for low electron density, the collision operator can be schematically written 



as 



Q[f] = [ [wik', k)/(k') - ^(k, k')f{k)]dk' (4) 

JR^ 

where the first term represents the gain and the second one the loss, and w{k, k') is 

the transition probability per unit time from a state k to a state k'. 

The main scattering mechanisms in a silicon semiconductor are the electron-phonon 

interactions, the interaction with impurities, the electron-electron scatterings and the 

interaction with stationary imperfections of the crystal, as vacancies. 

In this paper we shall consider electron-phonon interactions, and interaction with the 

impurities [IH|, i.e. 

W(k, k') = WaciK k') + WnpoiK k') + Wimp{k, k') . (5) 

The first term represents the acoustic phonon scattering in the elastic approximation 



wac{k,k') = Ko 6{e{k') - e{k)) 
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where Tl is the lattice temperature, the acoustic phonon deformation potential, p 
the silicon mass density and Vg the sound velocity of the longitudinal acoustic mode; 
the second term represents the non-polar optical scattering 



1 Tl'^pUJ 



where is the optical phonon energy and rig the phonon equilibrium distribution 
which, according to the Bose-Einstein statistics, is given by 

1 



exp{hhj/kBTL) - I 
The third and last term represents the ionized impurity scattering 



The values of the coupling constants and the other parameters used for silicon are 
given in Table I. 

The electron-electron interaction is taken into account in the framework of the 
mean field approximation through the Poisson equation. This is reasonable since we 
consider the case of low electron density and, therefore, we can neglect the short range 
collisions between electrons. 



3 Exponential tails 

In the following we shall draw our attention to the asymptotic behaviour of the steady 
state distribution function of the BTE 

U{k)= lim/(t,k) 

t^oo 

for large |k|, i.e. the high-energy tails. We expect that the behaviour of this solution 
is given by the so called 'stretched exponential', i.e. 

/oo(k) ^ exp [-re'^] , |k| ^ oo (6) 

where r and s are some positive constants. We introduce the functionals ^3] 

J'rAf)= [ fik)exp[re']dk (7) 

and study the values of s and r for which these functionals are positive and finite. 
We can give the following definition: 
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Definition. We say that the function / has an exponential tail of order s > 0, if 
the following supremum 

rt = sup{r > I J^rAf) < +00} (8) 

is positive and finite. 



The functionals ((Zj) can be represented by using the symmetric moments of the dis- 
tribution function, i.e. 

mp= I f{k)ePdk , peR+ . (9) 

In fact by expanding the exponential function in ((7j) into Taylor series we obtain 
(formally) : 



„ / 00 „ \ 00 



In order the expansion (fTUI) have a sense, we shall suppose that the moments of all 
orders are finite. Then the value r* can be interpreted as the radius of convergence 
of the series (jlOj) . and the order of the tail s is therefore the value for which the series 
has a positive and finite radius of convergence. For investigating the summability of 
the series fjlO|) we look for estimates of the sequence of moments {mp}, with p = sn, 
n = 0,1,2.. , and study the dependence of the estimate on s > 0. We shall be 
interested in the situation when the sequence of the coefficients satisfies 

^^CQ", n = 0,l,2... (11) 
nl 

where C, Q are positive constants. In fact, according to the Hadamard theorem 



hnsn 1 

lim r 



n— >oo y fi\ 

and from the estimate (|TT|l we have 

r:^i>0 . (12) 

To achieve the estimate (fTTjl . we shall study the moments equations obtained by 
integrating against the BTE, in the stationary and homogeneous regime : 

Qp + Gp = (13) 

where 

Qp = [ Q[f]e'dk (14) 



Gp = I E ■ VJe^dk (15) 
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We underline that in the stationary and homogeneous regime, the electric field E is a 
constant. Now we can formulate the main result of this paper: 

Theorem 1. Let /(k) be a nonnegative solution of the BTE ([T]) in the sta- 
tionary and homogeneous regime, that has finite moments of all orders. 
Then the supremum r* defined in (1121) is finite for some s> 0,and in the 
energy range huj < e < Eq, Weq > 0. 

To the best of our knowledge this is the first theorem which proves the existence 
of energy tails for the BTE in the semiconductor framework. This theorem proves 
also that the 'stretched exponential' (jH)), used in an heuristic way to model the EED 
tail, has a theoretical background. On the other hand, by using this technique, no 
information is given about the order of the tail s. This order will be determined 
by MC simulations in sec. 6. Our result is similar to that found in J3] for the 
BTE for granular media, where a shear fiow model is considered as forcing term. 
But the scattering mechanism for granular fiows is completely different respect to 
the semiconductor case, because binary inelastic collisions are considered between 
perfectly smooth hard sphere granular particles. The only analogy between the two 
cases, is the presence of an inelastic scattering mechanism in the collisional operator. 

4 Moment inequalities 



Let us introduce non dimensional variables 

7 k V £ ~ E Hlu 

k = — , V = — , e = — , E = — , u = — , a = aksTi 

ko Vo Bo Eq Eq 



where 



^/m*kBTL ksTi , „ „ kBTL^/m*k^rl 

ko = T ) ^0 = \/ — , £o = ksTi , Eq = 

n V qh 

consequently by omitting the" we can write 



k=2E{l + aEl v = ^^, . = . 

Now we change the variable k into spherical coordinates, i.e. 

k = {kx, ky, kz) — > {k sin 6* cos 0, k sin 6* sin 0, k cos 9) (16) 
and the density-of-state and the solid angle over the unit sphere are respectively: 

J = (1 + 2a£:) v^2e(l + aE), da = sin 6 d6 d(j) . 
The symmetric moments write: 



mp= / f{k)EPJdEda. (17) 
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Under suitable conditions on smoothness and decay for e large of the solutions /(k), 
the eq.(fT3j) writes: 



Jrs |E||k| e[l + 2ae) 

and from the definition of scalar product, and the eq. p7j) we have the following 
bounds : 

- 2j9|E|mp_i ^Gp^ 2j9|E|mp_i . (18) 
Now we manipulate the term Qp (I14|) : 

Qp= ! [w{V!,k)f{V!)-w{k,V!)f{k)\ePdV!dk= [ w{k,k')f{k)[{e'Y-eP]dk'dk 
and from the scattering rate (0) we have: 

Qf^ = [ 6{e'-e)f{k)[{er-{ey]dk'dk 

Q(npo) = Xi [ f{k)[ng6{e' ~e-u) + {ng + l)6{e' -6 + uj)][{e'y-{ey]dk'dk 
Qiimp) = K2 I 6{e' -e)f{k)[{e'Y -{eY]d\^dk 

The terms Q^^"* and Qp™''^-' vanish due to the presence of 6{e' — s). By changing 
variables into spherical ones, after some manipulation, we have: 



Qp = Qi""^ = Qi+^ + Q^-^ 

= 47T2-2n,K, I /(k) [{e + uof - e^] [1 + 2a{e + u)] 

X A/(e + uj)[l + a{e + 00)] Jde da (19) 



where 



1 



= 47r2t(n, + 1)7^1 / /(k) [{e - - e^] [1 + 2a{e - uo)] 
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x^{e-uj) [l + a{e-uo)]Q{e-uj)Jdeda. (20) 
Now we give an upper bound for Q^p!^\ By using the following inequalities 

yje + uj < e + 2uo 
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(where uj ~ 2.4369 in non dimensional units) and the definition p7|) . we obtain: 

+4cu(2a2cu2 + 3acu + l)mp + 2^'-^a2^*'+^m3 + 2^-^ujP{'ia + 8«^cj)m2 

+ (2tu)P(10a^cj2 + Qauj + l)mi + 2*^+2 (2a^w^ + Saw + l)mo} . (21) 

Now we give an upper bound for Q^p'\ By using the inequahties: 



\/e — UJ > a{e — oo) ,for oo < e < lo -\ — Va > 

as 

- (£ - iuY > e^~^, for e > Co- and p ^ — ~ 0.41 

we can write 

[{e - cu)^' - e^] < -£P~^a(£ - u f 

and Qp ■* can be estimated as 

^ -47r25(n, + l)iria / /(k)£P-i(£ - cu)^ [1 + 2a(£ - cu)] 

X a/1 + a{e — uj)Q{e — uj)Jde da. 
By using the further inequahties 

v/l + a(e-cj) > Pe, withe < ^ ,V/?>0 

1 + 2a{e - u) > Te, with e < \ ,Vr > 2a 

we obtain the final upper bound for Qp 

Ql'^ ^ -7r(ng + l)Kia(3r |2-imp+4 - 3cj2"^mp+3 + 3u}^2^mp+2 + ^^2imp+i| .(22) 

Now we are ready to obtain a recursive inequality for the moments. From eq. (jl3j) 
and the inequality (fTH|) . we get 

-2p\E\m^_^ <: Qp = Q; + Q; 

and using the upper bounds eqs.(j2IJ,(I22I), we obtain finally 

.1 ^ i?2lP-4lmp + 522lP-4lmp_i + 5i2lP-4lmp_2 + 5o2lP-4lmp_3 + T(p - 3)m 7 + 



y32b-4l2P-4^p-2^3 + 1/22lP-4l2P-4a;P-3m2 + Vi2lP-4l2P-3a;P-3^^ + 
^2lp-4l2P-icuP-2rno 
L Vp ^ 3 + ^ ~ 3.4104, < £ < £m 



where Em is an arbitrary fixed energy value, and the coefficients are positive constants 
given in the Appendix. 

Assuming some properties of the moments of lower order, we can use the recursive 
inequality to obtain information about the behaviour of the moments rrip, for p 
large. 



5 Proof of theorem 1 



We want to prove that the coefficients of the series (fTUjl satisfy the inequality (fTTjl . 
First of all we show that for every p ^ 1 there exist C and Q positive constants, 
such that 



nip ^ CQP . 

To obtain the desired estimate, we split our proof in two parts 

1. We fix po and we prove that the inequality (j^^ holds for 1 ^ p ^ po- 
We introduce the Jensen's inequality [T^ : 

Let be C[g{k)] a concave function and /(k) a density probability, then 



(24) 



f{k)C[g]dk <: C 



f{k)g{k)dk 



(25) 



For 1^ p ^ pq we choose: 

gik) = e'\ C[g]=g^o 

and from ()25|1 we obtain 



/(k) dk < J{k)sP''dk 

which can be written as: 

mp ^ (mp^)^ = CQP 
1 

with C=l and Q=max{l, (mpQ)po}. 



(26) 



2. We prove the inequality (j24|) for po < p < po + 1. 

We can estimate the moments /Wpo-i'^^po-ai'^po-S'^^po-- which appear in the 
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left-hand- side of by using the inequahty (^Bj). supposing p = Po ^ y ~^ h- 
The eq.(|^^ writes: 



where 



ii- 2? 



> 



Now eq.(j2Zj) can be written as 



PO 



(27) 



(28) 



where F = 2 max(ri,r2), and Q = max(2Q,4ci;) 

By using induction arguments, from eq. ()28|) we can conclude that the inequality 
(Plj) holds for po < p < po + 1. 

For every p ^ 1 we can cover the whole interval by using the previous arguments. A 
similar theorem is valid in the parabolic band approximation. 



6 Simulation results 

In the following we shall consider the bulk, where an homogeneous electric field is 
frozen into an homogeneous piece of silicon. The experimental evidence of the high- 
energy tails of the EED has been widely proved in literature by MC simulations. In 
the case of analytic parabolic band approximation (with acoustic and optical phonons 
scattering) this tail is a global maxwellian |Sl E] • So in principle we know, in this 
case, that the order of the tail in eq.® is s = 1. What is unknown is the radius of 
convergence r* of the series which is the inverse of the maxwellian temperature. 
This parameter will be determined by means of MC simulations. 
Several fitting formula are given in literature for the whole EED in the quasi parabolic 
band approximation. Here we shall check the EED tail by means of eq.® determin- 
ing, by using MC simulations, r* and s in the quasi parabolic band approximations. 
An histogram of the numerical steady solution is obtained by changing variable from 
k to {e,9,(f)) according to (fTH|) . discretizing the whole energy space in a system of 
concentric shells with increasing radius 

R 

Pn = nh, n = l,...,N, ^ = 
and by counting the number of particles which are in the corresponding shells 
/i = {# particles :£< pi} 

fn = {# particles : pn-i ^ e < pn} n = 2, N 
In+i = {# particles : i? ^ £:} 
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In order to obtain the parameters in ®, we rewrite this equation into the form: 



fn = cexp 



M 



,n ^ no 



and plot the pairs (a;„, i/n) 

Xn = lia{en - Sno) , 1/n = ln(ln - In /„) ,n = uq + 1, N 
Thus we expect that 

yn = SXn + b 

where 



b = —s InT; 



M 



(29) 



(30) 



is an almost linear plot for large values of the energy, whose slope is the exponent s, 
and the y- intercept is related to the parameter Tm (whose dimension is an energy). 
In figure 1 we show the simulation results obtained in the parabolic band approxi- 
mation (i.e. a=0 in eqs.®,©) with an electric field of 80 KV/cm. The solid curve 
shows the MC data while the dashed straight line represents eq.(|^ with s=l 
and Tm= 4907 K. Thus the asymptotic EED is the global maxwellian 



fo 



exp 



M 



In figure 2 we show the slope s as function of the electric field, both for parabolic 
and quasi parabolic band approximations. In the parabolic band approximation this 
value is approximately one (i.e. the EED has a global maxwellian tail), but in the 
quasi parabolic case the EED tail is given by the 'stretched exponential' The 
parameter s is a function of the electric field, and for small field values, approaches 
the value one, i.e. the tail is a global maxwellian. 

In figure 3 we show the parameter Tm as function of the electric field, both for 
parabolic and quasi parabolic band approximations. For moderate values of the 
electric field this parameter coincides with the lattice temperature (dashed straight 
line), but for higher values this parameter increases. This behaviour can be justified 
as follows: r* is the convergence radius of the power series ()10p. which is subject to 
the inequality ()12j). where Q is a function of the electric field E. 



7 Conclusions 

In this paper we have proved the existence of high-energy tails for the EED of the BTE 
(in the analytic band approximations, and scattering with acoustic and optical phonon 
and impurities), with a moment inequality technique, adapted from the granular flows 
|14j . We have also proved that the 'stretched exponential' ©, used in an heuristic 
way to model the EED tail, has a theoretical background. Our theoretical results are 
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confirmed by MC simulations for bulk silicon. In particular, in the parabolic band 
approximation, the tail is a global maxwellian whose temperature coincides with 
the lattice temperature for moderate electric fields whereas, for higher values, this 
temperature increases. In the quasi parabolic case the tail is a 'stretched exponential' 
(jni). In the future we shall try to extent the same technique to other scattering 
mechanisms, such as electron-electron. Another line of research will be to use other 
analytic techniques in order to obtain exact information about the order s of the tail. 
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Figure 1: Logarithmic plot of the EED obtained with MC simulation (solid curve) 
and with eq. ()29|l (solid dashed straight line) with s=l and Tm = 4907 K. The electric 
field is 80 KV/cm and the parabolic band approximation is used. 
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Field (V/cm) 

Figure 2: The parameter s as function of the electric field, for parabolic and quasi 
parabolic band approximations. The dashed straight line (s=l) indicates a global 
maxwellian tail. 
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Figure 3: The parameter Tm as function of the electric field, for parabolic and quasi 
parabolic band approximations. The dashed straight line indicates the lattice tem- 
perature 300 ° K . 
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Table I. Silicon parameters 



rUe 


electron rest mass 


9.109510-28 g 


m* 


effective mass 


0.32 me 


Tl 


lattice temperature 


300 


Po 


density 


2.33 g/cm^ 


Vs 


longitudinal sound speed 


9.18 10^ cm/sec 




acoustic phonon deformation potential 


9 eV 




optical pliouou energy 


63 nieV 


DtK 


optical phonon deformation potential 


11 10^ eV/cm 



Appendix 



n 2x 2^|E| 22 „ x 

°= t: — ,T = — ^ ,02 = — (3a + 8au!) 



-(lOa^cu + 9au + 1) , ^0 = uj{2a'uj^ + 3aa; + 1) 



2^ n 



2^ n„ 



aV(3 + 1 ar/? + 1 



lOa^a;^ + 9aa; + 1) , ^ = 



ar/3 + 1 



ar/3 rig + 1 



(2aa;2 + 2>auj + 1) 



X = 2q;2 + Mr/?a , and M e]0, 168]. 
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